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ABSTRACT 


1. Introduction 


Transportation plays an important role in the economic, production, and service systems, and it is a significant part 
of the Gross National Product (GNP) for each country. Therefore, researchers have improved routes, eliminated 
unnecessary trips, and built alternative shortcuts. The Vehicle Routing Problem (VRP) has tried to optimize the 
products distribution system to optimize related costs and other organization purposes [1]. The progressive 
urbanization development has made industries, especially support industries, the movement of people and 
commodities a problem whose complexity is constantly increasing. Each distributor wants to make the most of it, 
and on the other hand, they face problems such as traffic congestion, air pollution, wasted time on daily travel 
routes, increased fuel consumption, and vehicle depreciation to distribute commodities and gain liquidity [2]. 
Transportation has irreversible effects on the environment. For example, we can mention risks such as resource 
consumption, land use, toxic effects on ecosystems and humans, noise, and greenhouse gas emissions and 
pollutants. In addition, the emission of greenhouse gases and carbon dioxide is directly related to people's health 
and indirectly to the ozone layer degradation. Greenhouse gases emitted by transportation cause a significant 
amount of air pollution in different countries of the world [3]. Among the different transportation methods, 1/3 of 
greenhouse gases generated by transportation are for diesel vehicles for medium to heavy works [4]. The use of 
electric vehicles is one of the methods for reducing pollution. They have low emissions and play a considerable role 
in reducing emissions. Due to the fossil fuels limitation, the use of electric vehicles plays an important role in 
reducing fossil fuels consumption. Electric vehicles are also more efficient than diesel vehicles. 75% of energy in 
diesel vehicles is wasted as heat and friction, and only 25% of energy is converted to wheel drive force, while only 


20% of energy in electric vehicles is wasted. The number of parts used in electric vehicles is less than that of 
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conventional vehicles, and the heat generated is also low, so electrical vehicles have less maintenance cost than 


conventional vehicles. 


In some countries, governments provide facilities such as reducing taxes and duties on buying electric vehicles and 
giving them facilities. As another advantage, electric vehicles can be charged with electricity generated by various 
resources such as wind energy, solar energy, nuclear energy, water energy, and biofuels. The use of these resources 
reduces the dependence on oil and gasoline and, as a result, fuel imports are reduced, and it reduces fuel costs. The 
use of these devices is increasing in the transportation industry (public, commodity, etc.), and many studies have 
been done on the design and optimization of this type of vehicle. The main limitation in using electric vehicles is a 
low volume of battery energy compared to fossil-fuel vehicles, long-recharging time, and the limitation in refueling 
(charging) stations. The refueling in this vehicle is performed by unique equipment, and unlike other gasoline and 
diesel-fuel vehicles, fuel stations for them are limited. Hence, the limitation of charging stations must be considered 


for this vehicle. 


Another limitation is the dependence of energy consumption on the vehicle load amount that causes limitations in 
the real world. Due to their positive effect on air pollution reduction, the design and generation of refueling stations 
are increasing. The reduction of fuel consumption leads to a decrease in service cost and thus customer satisfaction 
[5]. The present paper focuses on the location-routing problem in electric vehicles. For that, we consider 
simultaneous multi-depots; the amount of energy consumption depends on the load in electric vehicles, 
simultaneous pick-up and delivery, and hard and soft time windows. For this purpose, we proposed a two-objective 
fuzzy-mathematical programming model for electric vehicles with a limitation in charge stations, the dependence 
of energy consumption to vehicle load, and simultaneous delivery and pick-up. We used the multi-objectives 
particle swarm meta-heuristic algorithms based on the Pareto archive and the NSGA-II algorithm to solve this 


model. 
2. Literature review 


Erdogan and Miller-Hooks, in 2012, investigated a transportation problem considering the fuel stations. Their 
proposed model was an integer programming model in which several virtual nodes were defined for each fuel 
station. That led to an increase in the problem dimensions and the solving solution. In the proposed model, vehicles 
were considered homogeneous, and the time window equations were used. To solve the problems, they proposed 
two new large-dimension heuristic algorithms. The first algorithm was based on Clarke and Wright savings 
algorithm, and the second was based on the distance density clustering algorithm. The computational results of the 
two algorithms were almost similar [5]. In 2014, Schneider et al. proposed a routing model considering electric 
vehicles, time windows, and fuel stations. In their model, virtual fuel stations were added. To solve, they applied a 
combined algorithm consisting of a Local Search (LS) algorithm and a Tabu Search (TS) algorithm. The numerical 
results obtained in the sample problems were reasonable [6]. In 2014, Philip et al. proposed the electric vehicle 
routing problem considering the technology of charging at the refueling stations. They also considered the half-full 
charge in the model. They used Simulated Annealing (SA) and local search algorithms to solve problems. In 


observations with 200 customers, local search algorithm, and more than 200 customers, the simulated annealing 
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algorithm had better performance [7]. In 2015, Goeke and Schneider proposed an electrical vehicle routing problem 
with a time window and considered a mixed fleet that included commercial electric vehicles and conventional 
vehicles. In this model, the fuel consumption rate was assumed to be a linear function of the traveled distance. Also, 
fuel consumption equations, speed, cargo slope, and weight were imposed in the limitations. A neighborhood 


search algorithm was used to solve them [8]. 


In 2016, Kog and Karaoglan introduced an integer linear programming model where there was no need to add 
virtual nodes for fuel stations. They modeled the green vehicle routing problem using the definition of auxiliary 
decision variables. They solved the proposed model using the simulated annealing meta-heuristic algorithm and 
branch and bound exact algorithm [9]. In 2016, Catay and Keskin removed the full charge limit, allowing the 
vehicle to make a full charge. This is more practical in the real world and reduces charging time. They proposed 
integer programming and used adaptive neighborhood search to solve. Their computational results showed that 
half-full recharging might improve routing decisions. Furthermore, the computational results showed that their 
proposed method is more effective in finding solutions compared to other similar research [10]. In 2016, Wen and 
Roberti considered the routing problem to visit each customer. The time window is included in the proposed model. 
The proposed integer programming model solved the problems with 46 customers reasonably. They proposed a 
three-step algorithm based on variable neighborhood search and dynamic programming to obtain quality solutions 
at the right time. Their computational results show that the proposed algorithm calculates small-size observations in 


0.1 sec and generates good solutions for more than 466 customers [11]. 


In 2016, Hiermann et al. proposed an integer linear programming model for routing problems with mixed transport 
fleets for electric vehicles. This model assumes that the capacity and cost of using vehicles vary depending on their 
type. In the model, they considered the time window and charging time limits of each vehicle. They used the branch 
and price algorithm, a heuristic combination method including large neighborhood search and local search, in 


addition to the notation process to solve the mode [12]. 


In 2016, Doppstadt et al. proposed a vehicle model with combined combustion and electric engines. They 
considered the transportation cost to each node in four states. In the first state, the vehicle is driven by the 
combustion engine, while the electric engine is charging. The vehicle is driven by combustion and electric engines 
in the second state. To solve, they used a heuristic 4-step method. Here an initial solution is generated first. Then 
they used local search and OPT-2 operators in the improvement phase [13]. Schiffer and Walther, in 2017, have 
studied the electric vehicle routing problem by considering the time window. They simultaneously studied the 
vehicle routing and the location of charging stations and considered charging different options in attention to the 
real world limitations. In their model, travel distance, vehicles number, charging stations number, and the total cost 
was minimized [14]. Camilo Paza et al., in 2018, studied the electric vehicles location-routing problem by 
considering multi-depots and time windows. They considered the availability of two energy supply technologies, 
including the conventional charging technology “Plug-In” and the battery replacement stations. They proposed 
three models. In the first model, conventional charging, in the second model, battery replacement stations, and in 


the third model, both technologies are considered [15]. Rohmer et al., in 2019, have studied the vehicle 
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* 
location-routing problem by considering time windows in the delivery system of the final perishable commodities 
[19]. Almohana et al. (2019) investigated the electronic vehicle routing problem with time windows in uncertain 
conditions. They presented a single-objective mathematical model and used meta-heuristic algorithms to solve the 
model [20]. Zheng et al. (2020) studied the vehicle routing problem in multi-warehouse multi-travel states with 
time windows and distribution time. They aimed to design a set of travels for the fleet of vehicles between different 


warehouses, to minimize the total travel time [21]. 


Ghobadi et al., in 2021, studied the electronic vehicles routing problem with fuzzy time windows and multiple 
depots by considering simultaneous pick-up and delivery. This paper presents an electric vehicle routing problem 
in a state of multi-warehouses with charging stations by considering an expected penalty for fuzzy time windows 
on delivery/delivery. Since this problem with fuzzy time windows and delivery/delivery limitations is an np-harp 
problem, three meta-heuristic algorithms, including simulated annealing, variable neighborhood search, and a 
combination, were used to solve the model [22]. Much attention has been paid to the electric vehicle 
location-routing problem in recent years. Whereas in previous research, Camilo Paza et al., in 2018, and Ghobadi et 
al., in 2021, have considered multi-depots in their research, and in other research, only one depot has been 
considered the model. On the other hand, attention has not been paid to the dependence of energy consumption on 
load in electric vehicles and simultaneous pick-up and delivery in previous research. Hence, to cover the existing 
research gap, the electric vehicle location-routing problem is investigated in the present study by simultaneously 
considering the multi-depots, the dependent energy consumption amount on the load in an electric vehicle, 


simultaneous pick-up, and delivery, and hard and soft time window. 


Table 1. A summary of previous research 
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With attention to Table 1, little research has considered multi-depots and the dependent charge consumption 
amount on the load in electric vehicles location-routing problem. Also, most of the presented models are 
single-purpose and definitive. While in this study, electric vehicles location-routing problem with multi-depots and 


ISSN: 2582-5267 www.ajbsr.net 


AJ B S R Asian Journal of Basic Science & Research 
JOURNAL Volume 4, Issue 2, Pages 01-31, April-June 2022 


the dependent charge consumption amount on the load in a fuzzy condition has been investigated, and a 
two-objective mathematical model has been designed for this problem. In addition, in previous research about the 
electric vehicles location-routing problem, the simultaneous pick-up and delivery of commodities have not been 
investigated, while this assumption has been taken into account in the present study. Another highlight difference 
between the present study and previous research is the solving method. In the present study, Gomez software has 
been used to solve small-size problems, and meta-heuristic multi-objective particle swarm algorithms based on 
Pareto archive and NSGA-II have been used to solve the large-size problems. The variable neighborhood search 


method improves the solutions in the proposed algorithm structure. 
3. Mathematical modeling 


This study investigated the electrical vehicles location-routing problem by considering multi-depots and time 
windows. This problem consists of multi-vehicles, several charging stations, and multi-depots which charging 
stations and depots are located. Here, the charge is considered fully, and the fuel consumption rate depends on the 
load that vehicles carry in each route. It means that the fuel consumption amount in a vehicle will be different at a 
given consumption rate on an unload route and the same route with a load. Hence, unlike previous research, we 
have considered the fuel consumption of electric vehicles to depend on the load carried. Considering this 
assumption can play an important role in choosing the optimal route and certainly has an effective role in visiting or 
not visiting fuel stations. In this model, simultaneous pick-up and delivery are assumed, and each point has two 
types of pick-up and delivery requests that are considered fuzzy. Also, service to demand points must be performed 


in a predetermined time window. 

The assumptions of the studied model are as follows: 

e Vehicles have a limited capacity. 

e The charge consumption by each vehicle depends on the vehicle load. 

e Some parameters are considered definite in this model, and some are indefinite. 


e The combined routing problem includes customers, multi-warehouses (depots), and a set of stations to charge 


the battery. 
e The location of depots and charging stations is considered. 
e Each charging station is allowed not to be visited at all or more than once. 
e Vehicles are needed to meet all customers. All customers must be met. 


e There are two types of time windows, hard and soft. Services that are outside the hard period and within the soft 


period, [Lb;, Ub;], (in terms of delays) are subject to a specified and fixed penalty cost per unit time. 
e Service outside the soft period is not allowed. 
e Pick-up and delivery are simultaneously considered. 


e The model is as two-objective. 


ISSN: 2582-5267 www.ajbsr.net 


AJ B S R Asian Journal of Basic Science & Research 
JOURNAL Volume 4, Issue 2, Pages 01-31, April-June 2022 


3.1. Model indices 

A set of all nodes, including customers, depots, and charging stations 
N.: A set of customers, c is customer indice. 

Ns: A set of potential nodes of charging stations, f is indice of charging station. 
Na: A set of potential nodes of depots, d is the indice of the depot. 

K: Number of vehicles, k is indice of vehicle. 

3.2. Model parameters 

dj;; The distance between nodes i and j. 

¢,: The variable cost of using a k™ vehicle per unit of distance. 

f.: Fixed cost of using a vehicle. 

Co: Cost per unit of charge. 

Po : Charge consumption rate per unit of distance. 

a : Charge consumption rate depending on the vehicle load per unit of distance. 
W,: Weight per unit of each delivered and received product. 

d .: The amount of fuzzy demand delivered to the c customer. 

Pp c: The amount of fuzzy demand received from the c" customer. 
cap ,: Load fuzzy capacity of the k" vehicle (in kilograms). 

Q : The battery fuzzy capacity for any vehicle. 

tj: Travel time between nodes i and j. 

M: A selected large number. 

PI ;: Fuzzy penalty amount for any delay. 

[e,];]: Hard time window for i™ node. 

[Lb,,Ub,]: Soft time window for i™ node. 

f,: The establishing cost for a charging station at point f. 


fy: The establishing cost for a depot at node d. 


Lb; ej lL; Ub; 
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3.3. Model variables 

Ze: If the charging station is located at node f, it will be equal to 1. Otherwise, it will be equal to 0. 
Za: If the depot is located at node d, it will be equal to 1; otherwise, it will be 0. 

Z: If k" vehicle is used, it will be equal to 1, and otherwise, it will be equal to 0. 


Xo: If k™ vehicle, which started its travel from depot d, travels between two nodes i and j, it is equal to 1; 


otherwise, it is equal to 0. 


Ya: The amount of demand delivered to the customer in node j, if the k" vehicle travels between nodes i and j (the 


vehicle that started its travel from depot d). 

zijn": The amount of demand received from the customer in node j if the k" vehicle travels between nodes i and j. 
a;: The time which a vehicle reaches the location of each node i. 

p;: The time which a vehicle leaves the location of each node 1. 

St;: The service time for each node 1. 

Yex: The charge amount that the k" vehicle has when meeting the c" customer. 

Yn: The charge amount that the k" vehicle has when it reaches the f" charging station. 

wij‘: The load amount that k™ vehicle carries between nodes i and j. 

YI: The delay in servicing to node i 


3.4, The main structure of the mathematical model 


Nf Na K K Na N 
min z1 = > fy XZ + >. fa X Za + >. fcxZ+ » > Cadyxh? + Co X dij (Po + awit xi") (1) 
f=l d=1 k=1 k=1d=1 i=1 j=1 
The objective Equation (1): The minimizing costs 
n 
minz2 = > L,Y; (2) 


N K 
> dr a1 vi EN, (3) 


K 
> xk vieN (4) 
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Limitation (4) shows that vehicles entering a node will leave that node, certainly. 


N N 
>, wit = > aft vd € Ng, Vk EK (5) 


Limitation (5) shows that vehicles that started their travel from depot d will certainly return to the same depot at the 


travel end. 
N 
> xt s1 WiEN,,dE€Ng, VWkEK (6) 


Limitation (6) shows that a charging station may be met or not. 


N Na K N Na K 

kd ,,kd kd, kd _ J H 
P2222 VE OD) 
j=ld=1k=1 j=ld=1k=1 
N Na K n Na K 

kd ,kd kd kd _ = q 
j=ld=1k=1 j=ld=1k=1 


The two limitations (7) and (8) ensure that the delivery/pick-up request for every i customer are met. 


n K 
Dt kay yy west =>" wh ka pid 5 “a wk vie N,d € Na (9) 
jHik=1 


j=1lk= 


Limitation (9) calculates the vehicle load amount. 
ka ViLJEN, dEN,kKEK (10) 


Limitation (10) ensures that the vehicle load does not exceed its capacity. 


dij (Do + &wis)xk* — Q(1 — xh") S Vin — Vix Vie Nc,jEN.dENgkEK (11) 
Vir — Vie S dij Go + awh) xi + Qa — xf vi € Nc,j € Nco,d € Ng, k EK (12) 
dij (Po + &wi')xi* — Q(1 — xi") < vin —Vjx Wie Ne, j € Np,d © Ng, k EK (13) 
Vir — Vie S dij Go + awh) xi* + QC — x4 Wi € Ne, j € Nr,d € Ng, k €K (14) 


Limitations (11) to (14) calculate the residual charge amount at each node. 
Vir = dia(Do + Ewik)xis* Vie Nod € Ng k EK (15) 
Limitation (15) ensures that the vehicle is sufficiently charged to return to the depot. 
Vin = Q Vie NaUNy, VREK (16) 


Limitation (16) ensures that the vehicle receives charge according to its capacity when meeting the charging 


station, and each vehicle leaves the depots on a full charge. 
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a; = (a; + St) -((1-xh4)M) +t VkEK, ifENdeNg (17) 


ay < (a; + Sti) +((1-xh4)M)+ tj, VREK, ifEN,dENg (18) 


Limitations (17) and (18) are related to the time window and feasibility of the schedule of each node, and if node j 


is not serviced after node i, these two limitations are disabled. 


P, =0 (19) 
Lb; < a; < Ub; VieN (20) 
p; > St; +a; vVieN (21) 
a; — 1; < Yl; VieN (22) 


Number (19) shows that the origin leaving time is at the moment 0. Limitations (20) and (21) show that the arrival, 
leaving, and service times for each node must be within a soft time window. Limitation (22) express the delay 


amount in service for any node. 
xi4 <Mzqg WREK, ijEN,deNg (23) 
Limitation (23) ensures a vehicle starts its travel from a depot to cross the edges if that depot is established. 


xK4 <Mzp WkEK, iorjeENp,d €Ng (24) 


Vir S$ Mz, VEEK, i EN, KEK (25) 
Limitations (24) and (25) ensure that a vehicle meets a charging station if that station is established. 
xi Za, Ze = {0,1}, Vir Zi, VEO, Wik, Gi >0 VijJEN,adeNg (26) 
Limitation (26) shows the values range of the model decision variables. 


According to expressions (1) and (2), the first objective function includes minimizing the costs of locating depots, 
charging stations, and using vehicles. Minimizing these costs means making up fewer depots and charging stations 
and using fewer vehicles. In this state, the service level or in-time services decreases and causes delay. The delay 
increases with a decrease in the cost of locating depots, charging stations, and using vehicles. Decreasing the first 


objective function increases the second objective function, which these two objective functions conflict each other. 
3.5. Defuzzification, the model 


Many methods have been proposed to solve fuzzy mathematical programming problems. The method of ranking 
proposed by Jimenez et al. in 2007 is used in this research. Jimenez proposed a method of ranking fuzzy numbers 
by comparing their waiting intervals [16]. If a triangular fuzzy number is as A={L,M,U} , it can be written as 


follows (Fig.1): 


f(x) = — L=x=M 
wA(xX)=4 1 X=M (27) 
ga(x) = M<X<U 
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Fig.1. Triangular fuzzy number 

To ensure that the inverse of functions f A (x) g nes exists, it is assumed f A (x) to be continuous and ascending 
and § 4 (x) continuous and descending. The expected range of a fuzzy number is defined as follows: 


EN(A) =[ER ES =| f° xdt, C0) - [ads ,00 | (28) 


By aggregating the components and changing the variable, we will have: 


e1(A)=[87,E7 1=| [ te (ada -f g;' (adar| (29) 


If the functions f A (x) BA (x) are linear and A a triangular fuzzy number, the expected range will be as follows: 
a | 1 (30) 
EI(A) = ns a +U)| 


Also, the expected value of the fuzzy number A is half the value of the expected range. 


BE’ #E 
EV (A) = G1) 


Furthermore, for a triangular fuzzy number, A it is as follows: 


EV(A) = oe (32) 


Definition 1. For both fuzzy numbers A and B , and membership degree for A greater than B is as follow: 


0 ifE; —E, <0 
Ly, (A, B) = aE ing tpt _ pg? eB] 
Bh, aE) =k, ) (33) 
1 ifE} -E? >0 
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EY Es | ana [E; »Ey ands ends 
So that [ 19 | and [ 1o49 | are the expected ranges for Aand B. When ,,, (A,B) =0.5, Aand B are the 


same. 


When zz, (A, B) =>a@,A with the least, a degree a is greater than or equal to B , it can be written as A 24 B 
Definition 2. Suppose x < R” is a vector; it is acceptable with degree aif main{ z2,, (Ax. B)} = @ it can be 


written as AX a, B. Equation (31) can be rewritten as follows: 


ES = Ef 
ES —-E’ +E? -EB - 


(34) 


or 
[d-Q ES + a@Ef |x = aE? +d-aQyE? 


Thus, according to the above definitions, we can convert the fuzzy model into a definite and accurate model. As 


follows: 


MinEV (C)x 


SL: (35) 
xe{xeR"|Ax>, B,x>0} 


Now, according to the above definitions and using the mentioned method, we convert the proposed fuzzy 


programming model into a definite model like it: 


Definite model of objective functions: 


Nf Na K 
minzl=) fpxzp+) faXzat) fox z 
f=i d=1 k=1 


K Na N N 420M + oY 
t+) >> > Cady ack casa 
k=1 d=1i=1 j=1 
ae +2a¢@+al 
-——> wis )xi@ (36) 


“PLE + 2PLM 4+ PLY 
min z2 = > else eel Le (37) 


y 2 
i=1 
Definite model of limitations (7) and (8) with B-section: 


N Na K 


~ di + e dM + av 

U A 
oD, xf yt) > 2 xityit = -pto +p 4S wien, (38) 
1k=1 


j=1d=1k=1 


iMez 
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> >. ait zit - > xf aykd — (4 — gb Pt ne +B Di ina viEN, (39) 
j=ld=1k=1 j=ld=1k=1 
Definite model of limitation (10) with B-section: 
capk + capt capi! + capl 
< [a- ) an ale OT i? ViLjJEN, dEN,kKEK (40) 
Definite model of limitations (11-14) with B-section: 
at’ +a™ at ee — 
a,j((1 — py Pot po o os i= 6 )—.-— +8 ipa — t= ~p)-—_.— 
+ ov ‘ed 
+p ——](1 - Xij ) S Vik — Vix ViENc,jEN,dENGKEK (41) 
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pe tor eee a — xk) Wie Ne,j EN d E Ng k EK (42) 
M 4 ,U en M4 qu . M 
aiy((1— py 4g ETP pS FO yg OPE ye (gp Pte 
Qe ang” 
+B ———](1 - xf 1) < Vip — Vix WiENc,j ENp,dE Ng KEK (43) 
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Definite model of limitation (15) with B-section: 


ps tps py + py a+ +a™ aM tal 
Yin 2 dij((L - B)->— + B >> +10 -B)—,— + 8 — Wixi Vie Nd © Nak 
EK (45 ) 
Definite model of limitation (16) with B-section: 
L4QM M4 Qu 
Vix =A - pie Q +p Cee ViENGgUN;, VKEK (46) 
4. The solving method 


The present study is focused on modeling and solving the “multi-depots location-routing problem” in electric 
vehicles by considering a simultaneous pick-up and delivery and time window and a multi-objective particle swarm 
algorithm based on the Pareto archive was proposed to solve the model. To continue, the proposed structure of the 
combined multi-objective PSO method was considered to optimize the two objective functions in the model 
simultaneously. The aim of designing the above method was to achieve more global optimal solutions (or more 
Pareto solutions). To evaluate this algorithm, its output was compared with the NSGA-II algorithm. PSO was first 


developed by a social psychologist named James Kennedy and an electronic engineer named Russell C. Eberhart in 
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1995. This development was based on previous experiences in modeling group behaviors observed in many birds’ 


species. The general structure designed for the combined multi-objective PSO method is in the following 


N initial and feasible responses be 
generated 


An improvement trend be applied 
to initial responses 


Checking the feasibility of 
improved responses 


Determinig the initial values of 
pg and pi 


Creating an initial Pareto archive 
collection 


Are the stop Reportting the Pareto 
conditions archive collection as 
met? output 


Responses be 
updated 


The improvement trend be applied on new responses. 
Checking the feasibility of the responses. 


Updating the values of 
pg and pi 


Updating the Pareto archive 
collection 


N responses of higher quality and with more scatter 
be selected as the next set of repetition responses 


Fig.2. Diagram of steps of multi-objectives particle swarm algorithm based on the Pareto archive 


4.1. How to display the response 
Many matrices were used to display each solution in this research, and each matrice corresponded to a variable 
defined in the model. For example, a one-dimensional matrice with N; (number of charging stations) was used to 


display the variable z;, a one-dimensional matrice with Ny (number of depots) was used to display the variable zg, 


and a one-dimensional matrice with N, (number of vehicles) was used to display the variable z,. 
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4.2. Generation of the initial solutions 


A random approach is used to generate the initial solutions in this paper. To generate the initial solutions, the Z;, Za, 
and Z, matrices were randomly generated, and then the rest of the solution matrices (model variables) were 
quantified in attention to the model constraints. If the population size is N, whenever a solution was generated as 
described, if not repeated, it will be added to the solution population. This trend will continue until the solutions 
number in the population reaches axN (a is a number greater than 1). The solutions generating method stops after 


aXxN repetition. 


On the other hand, the solutions number in each algorithm repetition equals N. Hence; we have to choose N 
solutions from aN available solutions as initial production orders. In this paper, the selection of the solutions 
initial population was based on a rapid method of sorting non-dominated solutions, described by Deb et al. in 2002 
[17]. In this method, aN solutions designed by the mentioned algorithm are sorted and leveled. Each level number 
indicates the quality of the solutions at that level. For example, the quality of the solution in level 1 is higher than 
the quality of the solution in level 2. Then, for the solutions in each level, a scale called the crowding distance was 
calculated in proportion to that level. This scale for the solutions of each level indicates the scatter of the solutions 
of that level.In this paper, a criterion called C, is defined for selecting the initial responses, which is obtained by the 
following Equation [18] 


622 (47) 
crowding _ dis 


The above criterion is calculated for each of the available responses. 
Rank: Indicates the level number where the response is located. 
Crowding dis: The crowding distance of each response is proportional to the rank of that response. 


After calculating the above criterion for all the solutions, the solutions were arranged ascending (based on the value 
of C,), and the first N solutions with fewer values of C, were selected as the initial algorithm solutions. The use of 
the C, criterion is based on selecting solutions with higher quality and more scatter as the initial population. The 
solutions generated are improved as much as possible using the improvement procedure. How the improvement 


procedure works will be described in the next section. 
4.3. Improvement procedure 


In this paper, an improvement procedure is designed that applies to the solutions (particles) and improves those 
solutions (particles) as much as possible. The designed improvement procedure in this research varies based on 
variable neighborhood search. In this procedure, three neighborhood search operators are designed combined with 


the VNS structure. This combination is as follows [18]: 
{for each input solution s 
K=1 


While the stopping criterion is met, do 
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S1=Apply NSS type k 
S=Acceptance method (S,S1) 


If s is improved then 


If k=4 then 
K=1 

Endif 

End while} 


As shown in the above structure, after applying the neighborhood structure to the solution, the acceptance 
procedure is applied to the resulting solution and the previous solution, and one of the two solutions is selected as 
the next repetition solution for VNS. The acceptance procedure determines and selects the dominant solution from 
two solutions using non-dominated relations. First neighborhood search structure: In this structure, the indices f, 
and f, are randomly generated in a uniform range [1,N¢], and the values of the cells corresponding to the indices 
generated in the matrice related to the variable Z; are exchanged with each other. Second neighborhood search 
structure: In this structure, the indices d, and d are randomly generated in a uniform range [1,Na], and the values of 
the cells corresponding to the indices generated in the matrice related to the variable Zy are exchanged with each 
other. Third neighborhood search structure: In this structure, the indices k; and kz are randomly generated in a 
uniform range [1,K], and the values of the cells corresponding to the indices generated in the matrice related to the 


variable Z;, are exchanged with each other. 


4.4. Particle update 

Here, genetic algorithm operators were used to update the particles [18]. How to update the particles is as follows: 
xi" = (x! — pl) +(x} =p Jas, (48) 

In the above Equation: 


t+1 


x'!_: i" particle in (t + 1)" repetition (generation) 
x! : i" particle in t" repeating. 
t 


p, : The best response which i" particle reaches it. 


p. : The best-found response. 


~S 


: A neighbor generated by a mutation operator for x; . 


L 


at 
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‘-” : This symbol shows the crossover operator. 
>+* : This is a selection symbol. 


To obtain the i” response in the (t+1)" repetition, five responses are generated: two responses are obtained from 
x! t x! P t 
crossover operator between~’ and P; | two are obtained from crossover operator on“ and “* 8 and one obtained 


t 
from mutation operator on ae Finally, the response with higher quality and more scatter among these five 


t+1 t 
responses is selected as *i In fact, Ps and i are used in this Eq. as a guide to reach next repetition responses. 


Crossover operator: 


The crossover operator designed in this algorithm is a single point crossover operator that is applied to the variables 
Z;, Za, and Z,. After two parents are input to the crossover operator, an index for each variable is randomly 
generated in a uniform range, and then two children are formed from the crossover of the two parents. Mutation 
Operator: The mutation operator used in the above Equation for particle updating is the same as variable 


neighborhood search, which is fully described in the previous section. 


4.5. Updating p;and p’, 


For each i" particle, among the neighborhoods found for this response, if there is a better neighborhood than p;, p; 
will be replaced by it. Otherwise, it remains unchanged. Of all the responses found so far, if the best response is 


better than pg, it will replace pg. Otherwise, it remains unchanged. 
4.6. Updating the Pareto archive 


The solution method used in the present research is based on the Pareto archive. A set called the Pareto archive is 
considered in the proposed algorithm, which holds the non-dominated responses produced by the algorithm. This 
set will be updated each time the algorithm is repeated. The method of updating is that the generated responses in it 
are repeated, and the responses existing in the Pareto archive are poured into a response pool and leveled together. 
Then among these responses, the responses existing in the first level (non-dominated responses) are selected and 


considered as the new Pareto archive. 
4.7. Select next-generation responses 


In each repeating, the algorithm needs a population of responses. To select the next repetition population in this 
research, the existing responses in its population are repeated, and the new responses generated by the algorithm are 
poured together in a response pool. After leveling and calculating the crowding distance criterion for each response 
according to the level of that response, using the Deb rule in 2002, N responses with the highest quality and the 


most scatter are selected as the repetition population after the algorithm. 
4.8. NSGA-II algorithm 


In this research, in order to solve the model, in addition to the multi-objective particle swarm algorithm, the 
NSGA-II algorithm based on the Pareto archive is also used. The diagram of the proposed structure of the NSGA-II 


algorithm is as follows. 
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N initial and feasible responses be 
generated 


An improvement trend be applied 
to initial responses 


Checking the feasibility of 
improved responses 


Mutation rate and crossover rate 
parameters be set 


Creating an initial Pareto archive 
collection 


Are the stop 


a > Reportting the Pareto 
conditions met? 


archive collection as 
output 


Responses be leveled and the crowding 
distance criterion be calculated 


Parents be selected and mutation and 
crossover operators be applied 


New and previous responses in a set be leveled and the 
crowding distance criterion be calculated 


N responses of higher quality and with more scatter be 
selected as the next set of repetition responses 


The improvement trend be applied to the selected responses 
and the Pareto archive be updated 


Fig.3. Diagram of NSGA-II algorithm steps 


How to generate the initial responses in this algorithm is random. N feasible responses will be randomly generated 
as the first repetition population of the algorithm. To apply the crossover operator, first, the parents must be 
selected for the crossover. In the present study, two higher quality responses are selected to select the parents using 
the non-dominated Equation and the double tournament method. The NSGA-II algorithm has two operators, 
crossover and mutation. The mutation and crossover operators used in this research are the same mutations and 


crossover operators used in particle updates in the particle swarm algorithm. This similarity is that the comparison 
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of the two algorithms is impartial. The particle swarm algorithm and NSGA-II algorithm can be compared under 


the same conditions. 
4.9. Dominant relation 


Due to the conflict between objective functions in the present research, using a multi-criteria optimization 
approach to solve the problem is more appropriate than single-criteria optimization methods. Hence, descriptions 
of the principles and foundations of multi-criteria optimization are provided in the following. The main difference 
between multi-criteria and single-criterion optimization problems is the existence of different and conflicting 
objective functions in multi-criteria problems. The existence of such goals makes it impossible to achieve the 
optimal solution (or solutions) through conventional algorithms for single-criteria optimization problems. To better 
understand the nature of multi-criteria optimization problems, the following principal definitions are provided (in 


all definitions, the minimization model below is considered with p decision variable and q objective function. 
sy ER qd] — Min y=flx)=(fi(),f20), fol) 
Definition 1: Domination Relation 


In multi-criteria optimization problems, the vector +4 overcomes the vector Q , if the following two conditions 


are met: 


7.0, )=75) 1D cd 
4G, 20 SE sie 


Therefore, the main goal in such problems is to find a set of points that dominate other points. 
4.10. Sorting and leveling responses 


The NSGA-II algorithm sets and levels the response in the population based on non-dominated relations. 
Relationships defined in section 9.4 were used to set and level the responses. First, all responses were compared 
using the mentioned relations, and the responses that no response can overcome were considered first-level 
responses. Then, this trend was repeated for the set of unspecified responses, and the next levels were specified. A 


response is of better quality when its level number is low. 


Therefore, levels with lower numbers were first used to select the responses. If it is possible to choose between two 
responses on the same level, the crowding distance criterion is used. The priority of a response for selection is high 


when the value of the crowding distance criterion for that response in a level is high. 
5. Computational Results 


First, to check the validity of the model and algorithm, the model of a small-size sample problem was solved by the 
Gomez software and multi-objective particle swarm algorithm, and then their results were compared. After 
validating the model and algorithm to examine the performance of MOPSO and NSGA-II algorithms, we 


performed them in the MATLAB software. The results of its performance in experimental problems generated by 
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them were compared based on comparative indices of quality, uniformity, scatter, and solving time. All 


calculations were performed using the core i7 7500U-12GB-1TB -R5 M335 4GB computer. 
5.1. Comparative indices 


There are many different indices for evaluating the quality and scatter of multi-objective meta-heuristic algorithms. 
For comparison in this paper, three indices, including quality, uniformity, and dispersion, are considered [18]. They 


are described in the following. 


Quality indice: This indice compares the quality of Pareto responses obtained by each method. This indice levels 
all Pareto responses obtained by artificial bee colony algorithm and genetic algorithm together and determines what 


percentage of level 1 belongs to each method. An algorithm has higher quality when this percentage is higher for it. 


Uniformity indice: This indice checks the distribution uniformity of Pareto responses obtained at the boundaries 
of the responses. This indice is defined as follows: 
em roi |Gragpn = d;| (49) 

(N _ 1) x Amean 
In the above Equation, d; represents the Euclidean distance between two adjacent non-dominated responses, and 


dinean 1S the average of d; values. 


Scatter indice: This indice is used to determine the amount of non-dominated responses found on the optimal 


boundary. The definition of scattering indice is as follows 


N . . 
D= >, _max(ll ~ yil)) (50) 


In the above Equation, ||xi - yill represents the Euclidean distance between two adjacent responses x} and y} on 


the optimal boundary. 
5.2. Experimental problems 


In this paper, several experimental problems are selected from the benchmark data of the electric vehicle routing 
problem. Since the benchmark sample problems do not fit the model presented in this research, some parameters 
have been adjusted based on existing data and previous research to check the model validity. In this paper, three 
small-size problems based on the C101 problem were designed from EVRPTW sample problems and solved using 
the Gomez software and the proposed algorithm. To design the problems from problem C101, several stations are 
considered depots and a number as stations, which number of depots, number of stations, and number of customers 
is different for each problem with small size. 


Table 2. Small-size problems 


Number of customers Number of stations Number of depots No 
5 3 3 1 
10 3 3 2 
15 3 3 3 
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After evaluating the model validity and the solution algorithm, other problems in the EVRPTW benchmark with 
100 customers were selected and solved by a multi-objective particle swarm algorithm and NSGA-II algorithm. 


Then, the performance of the two algorithms was compared based on the described indices. 
5.3. Setting the parameters 
The algorithm parameters are set as follows: 


In a multi-objective particle swarm algorithm, the population size is considered 500, the number of repetitions in 
the VNS algorithm, 10, and the number of repetitions of the algorithm, 500. In the genetic algorithm, the rates of 0.8 
and 0.1 are considered for the crossover and mutation, respectively, and the population size is equal to 500, and the 
algorithm repetition number is 500. To set the model parameters, the parameters in the EVRPTW benchmark are 
used, and the other parameters in the model are set as follows: To generate triangular numbers related to each fuzzy 
parameter (m1, m2, m3), first, m2 is selected, then a random number (r) is generated in the range (0,1). The m1 and 
m3 will be generated using the Equation m2*(1-r) and m2*(r+1), respectively. To quantify, the fuzzy parameter m2 
(according to the available data, if any) is determined, and the m1 and m3 values are determined using MATLAB 
software. For this reason, for setting these parameters, we only need to mention the value of m2. The following 


values are considered in the generation of sample problems. 


In all periods for the demand amount of customer pick-up and delivery, m2, based on the data in the EVRPTW 
benchmark, is the same and equal to the existing demand. Charging capacity, vehicle load capacity, the distance 
between nodes, and charge consumption rate are considered based on data in the EVRPTW benchmark. Also, 
among the mentioned parameters, load capacity, charge consumption rate, and charge capacity is fuzzy as 
described. Travel time between nodes is proportional to their distance and is considered uniform [1....20]. The soft 
and hard time window is determined based on the data of ready time values and due date values in the EVRPTW 
benchmark. The ready time is considered equal to the low limit of the soft time window, and the due date is the up 
limit of the hard time window. Product weight is considered equal to 5. The costs of establishing a depot and a 
charging station are generated in a uniform range [1000000....10000000], and the fixed and variable costs of using 
the vehicle are estimated at 2000000 and 100, respectively. Finally, the delay penalty is a fuzzy number from (m1, 
m2, m3) that m2 is generated in a uniform range [100....500]. The value of B for ranking the fuzzy numbers is 


considered 0.8. 

5.4. Model validation results 

To evaluate the model validity, first, the two-objective model was converted into a single-objective model using the 
LP-Metric method. Then, the single-objective model was solved in the Gomez software environment for small-size 
problems. The LP-Metric method was used in the present study. It is one of the most popular methods for 
multi-objective problems. We always try to minimize the deviations of the objective functions from their optimal 


value. In this method, the individual responses for optimizing each objective function were first calculated, then the 


following objective function was minimized. 


minz = [w*Gi@%) — AG )/AG DI) +10 —w)RO) -AeV/AC)] (51) 
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In the above Equation, f\(x*) is the optimal value obtained from solving the model considering the first objective 
function. The f, (x) is the value of the second objective function based on the optimal solution of solving the model 
with only the first objective function. The f2 (x*) is the optimal value obtained from solving model considering the 
second objective function. The f; (x) is the value of the first objective function based on the optimal solution of 
solving the model with only the second objective function. Also, w* is the weight of the objective function. The 
proposed model was solved by GAMS coding software and BARON solver. Here, the value of P was considered 1, 
and the weight of the objectives was considered the same and equal to 0.5. Small-size problems are solved by the 
MOPSO algorithm and Gomez software to optimize the LP-Metric method. The results of Gomez were evaluated 
based on the variables values and the constraints. The study results confirmed the feasibility and validity of the 


model. The results of solving the model by the proposed algorithm and Gomez software are shown in Table 3. 


Table 3. Solving results for sample problems of small size 


NSGA-II MOPSO GAMS 


Time z frolx*) | fiG@") | Time | z | fo") | fi@*) | Time | Z | f(x’) f(x") 


0.05 | 0.183 | 7243.2 | 124440.5 | 0.065 | 0.367 | 8031.9 | 116300.6 | 0.001 | 0.476 | 3380.5 478213.3 


0.07 | 0.145 | 8397.9 | 187342.1 | 0.07 | 0.269 | 8263.5 | 194291.1 | 0.004 | 0.236 | 3876.4 480561.7 


0.07 | 0.284 | 8149.5 | 272664.7 | 0.07 | 0.284 | 8149.5 | 272664.7 | 0.012 | 0.456 | 4069.1 532906.1 


According to Table 3 and based on the objective function values, the output solutions for each problem were 
non-dominated relative to each other. In some problems, the objective function values resulting from the genetic 


algorithm and the particle swarm algorithm were similar. 
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Fig.4. Comparison of objective functions values for problems of small size 
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On the other hand, comparing the objective function value of the LP-Metric method showed that in some cases, the 
Gomez software performed better than the two algorithms, whereas in other cases, the performance of the two 
solution algorithms was better. Considering that the results of Gomez are non-dominated over the particle swarm 
algorithm and genetics algorithm and the solution of the two algorithms are non-dominated relative to each other, 
the quality of results obtained in the three methods was at the same level. Hence, it can be concluded that the 


solution algorithm can solve the model efficiently and converge towards the optimal solution. 
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Fig.5. Comparison of performance time for problems of small size 
5.5. Performance results 


In this section, 30 problems were selected among the experimental problems in the EVRPTW benchmark and 
solved using the MOPSO and NSGA-II algorithms. Results obtained from informing the two algorithms according 


to the comparative indices are shown in the following Tables. 


Table 4-a. Results of solving EVRPTW sample problems 


MOPSO NSGA-II 
Number Number 
Prob . A , : : ; : , 
Quality | Spacing | Diversit of Pareto | Quality | Spacing | Diversit of Pareto 
y | sp g y CPU y |p g y CPU 
metric | metric metric archive metric | metric metric archive 
TIME TIME 
responses responses 
C101 94.3 1.12 2686.4 | 220.7 87 5.7 0. 61 1511.2 80.7 60 
C102 97.1 0.94 2385.6 | 230.1 67 2.9 0. 70 1741.2 135.4 59 
C103 72.9 0.71 2923.6 | 247.1 77 27.1 0.44 2441.3 119.4 78 
C104 100 0.85 2256.7 | 223.1 85 0 0.60 1437.9 147.9 88 
C105 88.6 0.76 2854.1 284.4 93 11.4 0.33 1607.5 141.5 94 
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C106 100 0. 92 2250.2 | 219.5 95 0 0.57 1349.4 | 160.4 62 
C107 | 77.6 0.69 2472.2 | 222.6 79 22.4 0.38 2035.1 | 159.1 80 
C108 | 85.9 0. 86 2780.7 | 217.1 63 14.1 0.54 2152.1 | 103.4 76 
C109 | 98.7 0. 80 2966.3 | 222.8 63 1.3 0.57 1708.5 | 139.9 57 
C201 100 0. 75 2127.4 | 243.6 67 0 0.54 1696.5 89.8 70 


Table 4-b. Results of solving EVRPTW sample problems 


MOPSO NSGA-II 
Number Number 
Prob Quality | Spacing | Diversity eon of Pareto | Quality | Spacing | Diversity eng of Pareto 
metric | metric metric aTiMi archive | metric | metric metric IME archive 

responses responses 
rl01 | 73.8 0.78 3145.3 | 231.1 91 26.2 0.58 1652.9 | 133.7 63 
r102 100 1.2 2960.9 | 292.3 67 0 0.45 2110.8 | 144.4 89 
r103 100 0. 79 2614.1 | 243.1 90 0 0.41 2077.5 | 120.1 69 
r104 | 84.1 0.68 2553.1 | 218.5 67 15.9 0.56 1701.1 | 117.9 78 
rl105 | 93.8 0. 93 2566.1 | 290.5 94 6.2 0. 61 1035.6 | 160.5 64 
1106 | 73.4 0. 81 2397.6 | 297.9 71 26.6 0. 73 1211.3 | 130.9 81 
rl07 | 82.1 1.1 2640.2 | 243.9 65 17.9 0.64 1282.9 | 131.8 68 
r108 100 0.90 2642.9 | 211.1 67 0 0.67 1914.2 | 155.9 83 
rl09 | 93.6 0.88 3011.2 | 225.8 82 6.4 0.52 1877.7 | 150.5 85 
r110 100 0. 94 2983.8 | 240.9 76 0 0. 56 1684.4 | 127.7 87 


Table 4-c. Results of solving EVRPTW sample problems 


MOPSO NSGA-II 

Number of Number 

Prob | Quality | Spacing | Diversity | CPU Pareto Quality | Spacing | Diversity | CPU | of Pareto 
metric metric metric TIME archive metric metric metric TIME archive 

responses responses 
rc101 89.3 0.52 2803.2 259.5 71 10.7 0.46 1449. 2 88. 4 75 
rc102 70.1 1.22 2484.3 226.2 90 29.9 0.34 1550.7 93.9 60 
rc103 95.3 0.67 3003.9 260.3 80 4.7 0.43 1587.3 158.7 66 
rc 104 100 1.06 2669.4 271.1 719 0 0.64 1435.1 72.9 94 
rc105 90.1 1.1 2450.8 222.2 94 9.9 0.49 1528.3 118.9 63 
rc106 92.5 1.1 3156.8 211.7 68 75 0.49 1740.4 86.7 90 
rc107 92.1 0.56 3081.1 229.7 87 7.9 0.44 1835.9 167.8 719 
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rc201 81.1 0.76 2690.2 231.9 87 18.9 0.62 1283.6 141.2 97 
rc202 89.3 0.64 2776.9 242.4 72 10.7 0.45 1577.5 120.4 60 
rc203 74.3 1.3 2734.4 250.8 80 25.6 0.66 1337.6 117.1 75 


According to the above Tables, in solving the electric vehicles location-routing problem by considering 


multi-depots, time windows, and simultaneous pick-up and delivery, the MOPSO algorithm in all states has a 


higher ability than the NSGA-II algorithm to generate solutions with higher quality. The MOPSO algorithm can 


generate solutions with higher scatter than the NSGA-II algorithm. In other words, the MOPSO algorithm has more 


ability for exploring and extracting the solution feasible region than the NSGA-II algorithm. The NSGA-II 


algorithm generates more uniform solutions than the MOPSO algorithm. In other words, the NSGA-II algorithm 


searches the solutions region with more uniformity. 
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In addition, a comparison of the performing times of the two algorithms to solve the experimental problems is 
shown in Fig. 9. It shows that the problem-solving time by the MOPSO algorithm is longer than the NSGA-II 


algorithm in all cases, which means that the MOPSO algorithm needs more time to solve these problems. 
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Fig.9. Comparison of performing time (sec) 
5.6. Statistical analysis of comparing two algorithms 


To solve the proposed model in this study, a multi-objective particle swarm algorithm and NSGA-II algorithm were 
used. Comparing the results of solving the sample problems by two algorithms was performed based on 
comparative quality indices, scatter, and uniformity. In this section, by assigning the appropriate assumptions, we 


studied the differences between the results of the two algorithms based on statistical analysis. 
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Assumption 1: There is a considerable difference between the quality indices of the MOPSO algorithm and the 
NSGA-II algorithm. 


Table 5. Test results of the first assumption 


Standard Sample 
Average error Average Sample 
deviation volume 
1.8 9.8 89.7 30 MOPSO 
1.8 9.8 10.3 30 NSGA-II 
Standard Sample 
Average error Average Sample 
deviation volume 
1.8 9.8 89.7 30 MOPSO 
1.8 9.8 10.3 30 NSGA-II 


As shown in Table 5, the average of the two groups is very different from each other. Also, the statistic value is 
outside the confidence range. Thus, the assumption HO is rejected, and the assumption H1 is accepted. This means 


a considerable difference between the quality indice of the MOPSO and the NSGA-II algorithms. 
Assumption 2: There is a considerable difference between the quality indice of MOPSO and NSGA-II algorithms. 


Table 6. Test results of the second assumption 


Average error Standard deviation | Average | Sample volume Sample 
0.036 0.199 0.87 30 MOPSO 
0.019 0.105 0.53 30 NSGA-II 


95% confidence range 
; Freedom Level of Mean 
degree | significance | difference Lower Upper 
MOPSO | 22.712 29 .000 .82800 .7534 .9026 
NSGA-IIE | 25.147 29 .000 48433 4449 5237 


As shown in Table 6, the average of the two groups is very different from each other. Also, the statistic value is 
outside the confidence range. Thus, the assumption HO is rejected, and the assumption H1 is accepted. This means 


a considerable difference between the uniformity indice of the MOPSO and the NSGA-II algorithms. 


Assumption 3: There is a considerable difference between the scatter indice of MOPSO and NSGA-II algorithms. 
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Table 7. Test results of the third assumption 


Standard Sample 
Average error Average Sample 
deviation volume 
50.2 274.8 2702.3 30 MOPSO 
56.9 S157 1651.8 30 NSGA-II 
Standard Sample 
Average error Average Sample 
deviation volume 
50.2 274.8 2702.3 30 MOPSO 
56.9 S117 1651.8 30 NSGA-II 


As shown in Table 7, the average of the two groups is very different from each other. Also, the statistic value is 
outside the confidence range. Thus, the assumption HO is rejected, and the assumption H1 is accepted. This means 


a considerable difference between the scatter indice of the MOPSO and the NSGA-II algorithms. 


Assumption 4: There is a considerable difference between the performing time of MOPSO and NSGA-II 


algorithms. 
Table 8. Test results of the fourth assumption 

Average error Standard deviation Average Sample volume Sample 
4.5 24.8 241.1 30 MOPSO 
4.9 26.5 127.2 30 NSGA-II 

Freedom Level of Mean 95% confidence range 

degree | significance | difference Lower Upper 
MOPSO _ | 53.128 29 .000 241.01333 231.7353 250.2914 
NSGA-II | 26.256 29 .000 127.18000 117.2734 137.0866 


As shown in Table 8, the average of the two groups is very different from each other, and the statistics amount is 
outside of confidence distance. Thus, the assumption Hp is rejected, and the assumption H1 is accepted. This means 


that there is a considerable difference between the performing time of the MOPSO and the NSGA-II algorithms. 


6. Conclusions 


The present study proposed a mathematical model for the location-routing of electric vehicles by considering hard 
and soft time windows and multi-depots. Here, the capacity of the vehicle was limited, and the charge consumption 


amount was dependent on the load of the vehicle. Vehicles were fully charged at charging stations. In this paper, a 
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fuzzy-mixed-integer programming model was first presented, and the Jimenez method was used to convert the 
fuzzy model to a definite model. The mathematical model had two objectives, (i) minimizing costs and (ii) 
minimizing penalties for late delivery of commodities or services to customers. After model designing, a combined 
multi-objective particle swarm algorithm and NSGA-II algorithm were proposed for solving the model. In order to 
evaluate the validity of the model and algorithm, the presented model was performed and solved in the Gomez 
precision software for small size problems, and the obtained results were compared to the solution algorithms 
results. Then, the problem-solving results by the proposed algorithm were compared to the solution results by the 


NSGA-II algorithm based on comparative indices of quality, scatter, and uniformity. 
In general, the results of solving the model showed that: 


e The model evaluation and solution algorithms confirmed the model validity and had a feasible area. Also, the 


solution algorithms can solve the model effectively and converge towards the optimal solution. 


e Comparative results showed that the MOPSO algorithm in all cases has a higher ability than the NSGA-II 
algorithm to generate solutions with higher quality and more scatter. Also, the NSGA-II algorithm has more 
uniform solutions than the MOPSO algorithm. In addition, the results showed that the performing time of the 
multi-objective MOPSO algorithm is longer than the NSGA-II algorithm. 


e There is a considerable difference between the quality indice of MOPSO and NSGA-II algorithms. 
e There is a considerable difference between the uniformity indice of MOPSO and NSGA-II algorithms. 
e There is a considerable difference between the scatter indice of MOPSO and NSGA-II algorithms. 


e There is a considerable difference between the performing time of MOPSO and NSGA-II algorithms. 
7. Recommendations for future research 

e Designing the model in the full-charge mode and half-charge mode. 

e Considering the vehicles with combined combustion and electric engines. 

e Considering probabilistic and fuzzy parameters together to apply uncertainty. 


e Using the other algorithms such as Whale Optimization Algorithm (WOA) or Imperialist Competitive 


Algorithm (ICA) to solve the model. 
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